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TRACKING RAPID INTRACELLULAR MOVEMENTS: 

A BAYESIAN RANDOM SET APPROACH 

By Vasileios Maroulas and Andreas Nebenfuhr^ 
University of Tennessee 

We focus on the biological problem of tracking organelles as they 
move through cells. In the past, most intracellular movements were 
recorded manually, however, the results are too incomplete to capture 
the full complexity of organelle motions. An automated tracking al¬ 
gorithm promises to provide a complete analysis of noisy microscopy 
data. In this paper, we adopt statistical techniques from a Bayesian 
random set point of view. Instead of considering each individual or¬ 
ganelle, we examine a random set whose members are the organelle 
states and we establish a Bayesian filtering algorithm involving such 
set states. The propagated multi-object densities are approximated 
using a Gaussian mixture scheme. Our algorithm is applied to syn¬ 
thetic and experimental data. 


1. Introduction. Most plant cells display a striking phenomenon called 
“cytoplasmic streaming,” a process that has been recognized since the late 
18th century by Corti (1774). During cytoplasmic streaming, most subcellu- 
lar organelles move rapidly through the cell, resulting in constant mixing of 
the soluble components of the cytoplasm. The function of these movements is 
not known, although a potential role in better distribution of metabolites has 
been proposed in Shimmen and Yokota (1994). The movements are driven 
by myosin motor proteins [Shimmen (2007)] and appear to be necessary for 
normal growth of plant cells and ultimately the whole plant [Peremyslov 
et al. (2008), Ojangu et al. (2012)]. The molecular mechanisms that con¬ 
nect the intracellular movements with cell growth are not known [Madison 
and Nebenfiihr (2013)]. Better understanding of these cellular processes re¬ 
quires the targeted manipulation of the movements followed by quantitative 
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assessment of the resulting changes at the subcellular, cellular and whole- 
plant levels. Recent results have identified additional regulatory mechanisms 
that influence intracellular movements, although the precise nature of these 
mechanisms is still unknown [Vick and Nebenfiihr (2012)]. This is due, at 
least in part, to the astounding complexity of these movements and the 
technical difficulty of describing them accurately [Nebenfiihr et al. (1999), 
Hamada et al. (2012)]. 

Recent advances in molecular biology and fluorescence microscopy imag¬ 
ing have made possible the detailed observation of these intracellular dynam¬ 
ics and the acquisition of large multidimensional image data sets [Danuser 
(2011)]. Paredez, Somerville and Ehrhardt (2006) noted that these time- 
lapse observations reveal a large number of nearly identical particles that 
move with high velocities in close proximity to each other. Combined with 
the saltatory, or stop-and-go, nature of their motions, these features make 
automated tracking of these movements an extremely difficult task as dis¬ 
cussed in Nebenfiihr et al. (1999). As a result, most previous analyses have 
relied on manual tracking of a few individual particles, for example, 
Nebenfiihr et al. (1999), Gutierrez et al. (2009), Hamada et al. (2012), Logan 
and Leaver (2000), Collings et al. (2002). A full understanding of the obser¬ 
vations, however, requires accurate tracking of a large number of bright spots 
in noisy image sequences, which can be accomplished only by an automated 
algorithm that is able to analyze the data completely [Danuser (2011)]. This 
complete analysis will require reliable identification of organelle positions 
(coordinates) from the bright spots in fluorescent microscope images taken 
at different times and the correct linking of these positions into continuous 
movement trajectories over all time points. One benefit of such an algorithm 
could be the emergence of recurring patterns such as the recent discovery, 
based on manual tracking, that organelles preferentially pause their motions 
at microtubules [Hamada et al. (2012)]. Thus, it seems likely that a compre¬ 
hensive and accurate tracking algorithm will unearth additional regulatory 
events that in turn can be studied experimentally. Moreover, from a statisti¬ 
cal point of view, an automated tracking algorithm will reduce the bias since 
manual tracking depends solely on experts’ decision of linking the positions 
of bright spots at subsequent time points. 

Mathematical and statistical models that require knowledge from statis¬ 
tics, probability, scientific computing and statistical mechanics have been 
developed for reliably tracking multiple objects in space. There are a great 
number of studies addressing the problem of tracking multiple targets in 
various settings. A partial list of such works is Doucet, de Freitas and Gor¬ 
don (2001), Liu (2008), Goodman, Mahler and Nguyen (1997), Gilks and 
Berzuini (2001), Fortmann, Bar-Shalom and Scheffe (1983), Bar-Shalom and 
Blair (2000), Blackman and Popoli (1999), Liu and Ghen (1998), Maroulas 
and Stinis (2012), Vo, Vo and Gantoni (2007), Mahler (2007, 2003), Mahler 
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and Maroulas (2013). However, only a small number of multi-object mod¬ 
els have been considered for specific microscopy image data, for example, 
Smal (2009), Smal, Niessen and Meijering (2006), Sbalzarini and Koumout- 
sakos (2005), Jaqaman et al. (2008). Movement of subcellular particles in 
living cells poses a highly complex problem for automated tracking algo¬ 
rithms. Even at high magnification, the true position of a particle within a 
cell can be only measured to within 50-200 nm due to limitations in optical 
resolution, and given the inevitable image noise, it is likely that some or¬ 
ganelles are not detected. Moreover, not only can individual organelles move 
independently, they also can change their behavior rapidly, their paths are 
not static, and organelles in close proximity can display strikingly different 
behaviors [Collings et al. (2002), Nebenfiihr et al. (1999)]. Commercial au¬ 
tomated tracking algorithms such as Perkin-Elmer’s “Volocity” were some¬ 
times used to gain insights into overall movement patterns or derive average 
movement velocities; for example, see Peremyslov et al. (2008), Avisar et al. 
(2008). However, these algorithms often introduced mis-assignments in the 
tracks [e.g.. Figure 3A in Avisar et al. (2008)] and, therefore, cannot be used 
to obtain an accurate global view of organelle motility. 

In general, from a statistical point of view, tracking of multiple objects is 
an inherently difficult problem and consists of computing the best estimate 
of the objects’ trajectories based on noisy observations. The estimates are 
propagated by a posterior distribution which considers organelles’ dynamics 
and combines them with data. The greater the number of objects that are 
being tracked, the more complicated the tracking algorithm becomes. There 
are several techniques, for example, Kalman filters and their derivatives, 
particle filters, for addressing this problem statistically. The reader may 
refer to Gordon, Salmond and Smith (1993), Liu (2008) and the references 
therein. 

A popular approach to tracking is particle filtering. Smal et al. (2008) 
introduced a particle filtering algorithm for the tracking problem using mi¬ 
crotubule dynamics, which overall follow a priori known and fairly straight 
paths and can therefore be conveniently modeled. In general, the parti¬ 
cle filter approach is an importance sampling method which approximates 
the posterior distribution by a discrete set of weighted samples (particles). 
However, it is often found in practice that most samples’ contribution to the 
posterior distribution will be negligible. Therefore, carrying them along does 
not contribute significantly to finding an estimate. Hence, one may resam¬ 
ple the particles to create more copies of samples with significant weights 
[Gordon, Salmond and Smith (1993)]. However, even with the resampling 
step, the particle filter might still need a large number of samples in order 
to approximate accurately the target distribution. Typically, a few samples 
dominate the weight distribution, while the rest of the samples are in statis¬ 
tically insignificant regions [Snyder et al. (2008)]. Thus, some studies [see. 
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e.g., Gilks and Berzuini (2001), Maroulas and Stinis (2012), Weare (2009), 
Kang and Maroulas (2013)] have used an additional Markov Chain Monte 
Carlo step which helps to move more samples into statistically significant 
regions and thus to improve the diversity of samples. This extra step can 
improve estimates for multi-target tracking scenarios [Maroulas and Stinis 
(2012), Kang, Maroulas and Schizas (2014)], but at the price of adding an 
additional layer of complexity. 

In this manuscript, we attempt to avoid the technical algorithmic steps 
which depend on the specific nature of different applications. Instead, we 
create an automated statistical tracking algorithm for independently evolv¬ 
ing intracellular movements by considering a pertinent multi-object statis¬ 
tical framework. This framework adopts a Bayesian random set hltering 
technique. The key innovation in our approach is to conceptually view the 
evolving collection of organelles as a single set-valued state and the collec¬ 
tion of the experimental measurements as a single set-valued observation. 
A set-valued state contains not only the position of existing organelles but 
also the states of new biological entities which enter the tracking domain. 
Using Random Finite Set (RFS) theory and modeling the collection of or¬ 
ganelles and their corresponding experimental measurements as sets result 
in generalizing single-object filtering to a rigorous formulation of Bayesian 
multi-object filtering. Multi-object hltering, similar to the single-object case, 
consists of two stages, the prediction stage using modeled or experimen¬ 
tally derived dynamics, and the update stage using the observed data. Both 
these steps involve multi-object distributions which lead to the multi-object 
Bayesian hltering posterior distribution, 

f{X\Zi, f{Zt\X)fiX\Zi ,..., Zt-i), 

where X,Zi,... ,Zt are appropriate random sets, formally dehned in Sec¬ 
tion 2. 

The general multi-object Bayes hltering distribution, f{X\Zi ,..., Zt), is, 
however, computationally intractable in most applications and thus it needs 
to be approximated. In this paper, we consider a Gaussian mixture Cardinal- 
ized Probability Hypothesis Density (CPHD) approximation. The CPHD, 
hrst introduced by Mahler (2007), propagates two estimates, the cardinality 
distribution of a random set which yields an estimate of the number of ob¬ 
jects per time step, and the intensity of a random hnite set or otherwise the 
so-called probability hypothesis density (PHD) [Mahler (2003)]. The PHD 
is similar to the hrst-moment density or intensity density in point process 
theory; for example, see Daley and Vere-Jones (1988). The PHD hrst mon¬ 
itors multiple objects as clusters, and then attempts to resolve individual 
objects only as the quality and quantity of data permits. One could also es¬ 
timate the number of objects at a given time step using the PHD, however. 
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such an estimate is unstable when the experimental scene is highly dynamic, 
that is, with rapid entry and exit of organelles from the region of interest. 
A Gaussian mixture approximation of the CPHD was introduced by Vo, Vo 
and Cantoni (2007) whose algorithmic complexity was of the order 
where m is the number of data points (acquired positions of organelles) and 
n the true number of objects of interest. However, the cubic dependency on 
the number of data points is disadvantageous for our biological framework 
due to their large number. 

In our manuscript, we consider a Gaussian mixture GPHD based on the 
experimental fact that data are generated only when organelles are present 
in the tracking domain. A false alarm is generated in signal detection when 
a nontarget event exceeds the detection threshold. Our experiments did not 
suffer from any false alarm, and thus a pertinent approximation of the GPHD 
is established in Propositions 2.1 and 2.2. The associated algorithmic im¬ 
plementation cost reduces to the order of 0{mn), that is, the cost is linear 
with respect to the number of data and organelles. In brief, Proposition 2.1 
propagates the predicted cardinality and the predicted intensity estimate 
(PHD) of a random finite set which follows a Gaussian mixture density. Tak¬ 
ing into consideration a new random set of data (positions of organelles). 
Proposition 2.2 updates the two predictions by considering a Bayesian set 
formulation. The posterior PHD follows an appropriate Gaussian mixture 
whose components are derived with the aid of Proposition 2.2. 

A similar algorithm was analyzed in Mahler and Maroulas (2013) for the 
special case of monitoring two hxed objects that spawn several objects along 
their ballistic trajectories. These secondary objects fall under gravity, and 
thus they are not of tracking interest. Precisely, a distance criterion was 
computed to distinguish the two primary objects from the spawned ones. 
When this distance exceeded a certain threshold, the corresponding objects 
were declared debris and they were discarded. This assumption cannot be 
incorporated herein. Thus, in our framework, we relax this condition and, 
moreover, we incorporate several experimental biophysical features to un¬ 
derstand the unknown dynamics of organelles. For instance, based on the 
organelles’ acceleration data analysis (see Section 3), we discover that the 
acceleration follows a normal distribution with mean-zero. Assuming that 
the mass of the observed organelles did not change significantly between 
individual images (a valid assumption), we are able to deduce interesting 
results about the developed biomechanics within a cell. 

Section 2 focuses on the methodology that was followed to establish an 
automated tracking algorithm for organelle movement data. Definitions of 
the Gardinalized Probability Hypothesis Density (GPHD) and approxima¬ 
tion schemes are also presented. Section 3 describes the implementation of 
an appropriate version of the Gaussian mixture GPHD filter suited for the 
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biological data (synthetic and experimental). Section 3.2 describes the bio¬ 
physical conditions under which the experimental data were collected and 
the process of manual tracking. Finally, our results are summarized in Sec¬ 
tion 4 and a discussion for future research and developments is offered. 


2. Random finite sets and approximations. We motivate this section by 
considering first the problem of tracking only one object. Suppose that an 
organelle, whose state is x' at time t, moves following the dynamics below, 

( 2 . 1 ) xt+\=(f>t{x' ,ut), 

where ut is a randomly distributed noise and (pt ■ 1^^ x —>■ is a family 

of nonlinear, nonsingular functions. Let zi-t = {zi,Z 2 -, ■ ■ ■, zt] denote the data 
history up to time t and let ft\t{x'\zi-t) represent the posterior probability 
density function (p.d.f.) at a given time t. Furthermore, consider the poste¬ 
rior predictive p.d.f., ft+i\t{x\zi-t), which merely yields the probability that 
an organelle will move to state x at time t + 1 given the available data zi-j- 
Using the Chapman-Kohnogorov equation, the posterior predictive distri¬ 
bution is given by 

( 2 - 2 ) ft+i\tix\zi-,t) = j ft+i\tix\x')ft\tix'\zi:t)dx', 

where ft+i\tix\x') is the Markov transition density associated with the dy¬ 
namics expressed of equation (2.1). At given time t-|- 1, a new microscopy 
observation is collected, zt+i € Typically, the dimension of organelle 
states, N, and the dimension of data, M, are not identical, N ^ M. For ex¬ 
ample, the state of organelles involves their position on the xy-plane and the 
corresponding velocities, that is, = 4, whereas only the positions (M = 2) 
are available from the experimental data. The prediction (2.2) needs to be 
updated using the datum zt+i ■ The collected measurement is a function of 
the true organelle’s state perturbed by noise, that is. 


(2.3) 


Zt+i=r]t+i{x,Ct+i), 


where P^t+i is a randomly distributed noise, independent from vt, and the 
function rjt+i : x —>■ is a family of nonsingular, nonlinear trans¬ 

formations. Based on the Bayesian rule, the posterior p.d.f. at a given time 
t -|- 1 is given by 


(2.4) 


ft+l\t+l{x\zi:t+l) 


ft+i{zt+i\x)ft+i\t{x\zi-,t) 

f ft+i{zt+i\x)ft+i\t{x\zi-,t)dx’ 


where ft+i{z\x) is the likelihood function associated with (2.3) and the pos¬ 
terior predictive distribution, ft+\\t{x\zi:t), is defined in (2.2). 
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Remark 2.1. The widely-known Kalman filter is a special case of the 
Bayesian filtering formulation given in equation (2.4). Indeed, if one con¬ 
sidered that were linear and vt,'Wt were normally distributed, then 

equations (2.2) and (2.4) would enjoy a closed-form solution which would 
be the same as in the Kalman filter. 

On the other hand, our focus is on tracking multiple objects which move 
simultaneously. Motivated by the single object tracking framework described 
in equations (2.2) and (2.4), we consider a statistical framework which allows 
us to generalize the prediction equation (2.2) and the corresponding update 
equation (2.4), both suitable for tracking one object to pertinent equations 
for tracking one set of objects. We view for the first time in this biological 
problem the evolving collection of the organelles as a single set-valued state, 
Xt = {xl,x^,... ,xf*} € T'(M'^), where nt represents the number of objects 
at time t, and T'(M^) is the collection of all finite subsets of Similarly, 
the collection of experimental microscopy measurements at time t is viewed 
as a single set-valued observation, Zt = {z}, z^,..., ^ where 

nit is the number of generated measurements at time t. Based on equation 
(2.3), each member zl £ Zt+i is a noisy perturbation of the true state x of 
an organelle j at time t, where i is not necessarily equal to j. 

Furthermore, the randomness in this multi-object framework is repre¬ 
sented by modeling multi-object states, A^, and multi-object measurements, 
Xlt, as random finite sets (RFS) on the single-object state and observation 
spaces, and respectively. The corresponding multi-object dynamics 
and observations are described below. 

Given a realization, Xt, of the RFS, Aj, at time t, the multi-object state 
at time t -|- 1 is modeled by the RFS, 

(2-5) A^+i = I 5t_|_i|f(x)| U Rt+i, 

''3;GAt J 

where is the RFS representing the objects which survive with proba¬ 

bility Ps,t+i\t{x), from the previous time t, and Bt is the RFS which repre¬ 
sents the objects which enter the scene at time t-\-l (“newborn” organelles). 
Hence, the RFS, A^+i, includes all information of set dynamics, such as the 
number of objects that vary over time and an individual organelle’s motion 
[see equation (2.1)] and birth/death. Now, given a realization Xt+i of Aj+i 
at time t-\-1, the multi-object measurements are modeled via the following 
RFS, 


Mt+i — 0t+i(a;), 

xeXt 


( 2 . 6 ) 
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where 0t+i(x) is the RFS of measurements generated by the object x ^Xt- 
The RFS Mt+i encapsulates all characteristics of the measurements from 
the microscopy image, for example, measurement noise. 

Next, let denote the multi-object posterior density at a given 

time step t conditioned on the observation sets, Zi-t = {Zi, Z 2 -, • • •, Zt}. The 
multi-object Bayes filter propagates the multi-object filtering distribution 
via the following recursion: 


(2.7) 

( 2 . 8 ) 


ft+i\t{X\Z^:t) = I /t+i|i(X|X')/i+i|t(X'|Zi,i)5X', 


ft+i\t+i{X\Z^.,t+i) 


/t+i(Zt+i|X)/,+i|i(X|Zi,0 

f ft+iiZt+i\X)ft+,^t{X\Z,,t)6X^ 


where f 6X is the set integral [see, e.g., Goodman, Mahler and Nguyen 
(1997), Definition 10], ft+i\t{X\X') is the multi-object transition density 
associated with the dynamics given in equation (2.5), and ft+i{Zt+i\X) 
is the multi-object likelihood obtained by equation (2.6). One may show 
that densities and likelihoods expressed in equations (2.7) and (2.8) are well 
defined using techniques of finite set statistics (FISST) and extending the 
concept of the Radon-Nikodym derivative [Goodman, Mahler and Nguyen 
(1997), Ghapter II.5]. 


Remark 2.2. One may compare the analogy between equations (2.7), 
(2.8) and equations (2.2), (2.4), respectively. Therefore, our statistical frame¬ 
work generalizes the problem from tracking a single object to tracking a 
single set. 


However, the multi-object filter described in equations (2.7) and (2.8) is 
intractable in most applications and the Gardinalized Probability Hypoth¬ 
esis Density (CPHD) approximation is considered. The GPHD produces 
estimates on the number of organelles and their states. A formal definition 
is below. 


Definition 2.1. The GPHD filter recursively propagates the posterior 
cardinality distribution p^^^{n\Zl■t) on object-number n and the intensity 
function or Probability Hypothesis Density (PHD) Df^i{x\Zi-t)■ Given any 
region S C , the expected number of objects in S is derived by the in¬ 
tegral Jg Dm{x\Zi:t) dx. If S = M'^, then = f Dfn(xlZi:t) dx is the total 
expected number of objects in the scene. 

The GPHD filter produces stable (low-variance) estimates of object num¬ 
ber, as well as better estimates of the states of individual objects [Mahler 
(2007), Vo, Vo and Cantoni (2007)]. This gain in performance is achieved 
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with increased computational cost. For instance, Vo, Vo and Cantoni (2007) 
implemented a Gaussian mixture CPHD whose algorithmic cost was of the 
order 0{m?n), where m is the number of data points and n the number 
of objects of interest. However, the number of data points is large and the 
number of organelles is a priori unknown and varies in time. Therefore, 
the alternative Gaussian mixture implementation of Mahler and Maroulas 
(2013) is considered herein which decreases the computational cost to the 
order of 0{mn). In fact, our technique is based on the experimental obser¬ 
vation that all data are produced by the organelles and no false alarms exist. 
If false alarms were collected, for instance, due to human intervention, then 
equations (2.5) and (2.6) would need to be suitably formulated. 

Before proceeding with the dynamics and Bayesian formulations as ex¬ 
pressed in Propositions 2.1 and 2.2, respectively, we list the assumptions on 
which our Gaussian mixture approach to the CPHD is based. 

Assumption 2.1. Consider a realization Xt = {x\,xl,... of the 
RFS, At, and the associated data collection Zt = {z}, ,..., z^*}. The state 

of each organelle xl € Xt,i = 1,... ,nt is normally distributed given by 

(2.9) xt\xt-i ~ N{x;Ft-iXt-i,Qt-i), 

where Ft-i is the state transition matrix and Qt-i is the process noise co- 
variance. Similarly, each observation zlj = l,...,mt,j^i, is normally dis¬ 
tributed according to 

(2.10) ztlxt--^ N{z-,HtXt,Rt), 

where Fit is the observation matrix and Rt is the observation noise covari¬ 
ance. 

Assumption 2.2. The survival probability, Ps,t+i|t(x), of an organelle 
with state x at time t to be present at time t -|- 1 is state independent, 
that is, Ps,t+i\tix) =PS- The detection probability, PD\t+iix)j to collect an 
observation associated with an organelle whose state is x at a given time t, 
is state independent, that is, PD\t+iix) =Pd- 

Assumption 2.3. The intensity measure of the birth RFS which en¬ 
compasses the dynamics of newborn organelles is a Gaussian mixture of the 
form 

( 2 . 11 ) bt{x) = Y^w^^lN{x-,p'^l,Plt), 

i=l 

where the weights, means and covariances of the mixture 

birth intensity and Jh^t is the number of Gaussian components associated 
with the newborn organelles at a given time t. 
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Remark 2.3. Assumptions 2.1-2.3 are crucial for establishing a closed 
form for the multi-object densities dehned in equations (2.7) and (2.8). How¬ 
ever, if the linearity of Assumption 2.1 is violated, then one could consider 
implementing a CPHD hlter introduced by Vo, Vo and Cantoni (2007), which 
employs a pertinent approximation of the nonlinearities. However, Assump¬ 
tions 2.1-2.3 are satished using our experimental data. Further discussion 
of this topic is delegated to Section 3. 

The propositions below involve the main equations of the Gaussian mix¬ 
ture implementations of the CPHD hlter without considering any false alarms 
For presentation’s sake, the time index is suppressed in the cardinality of the 
state sets and measurement sets in the propositions below, that is, nt = n 
and mt+i = m. The reader should refer to Mahler and Maroulas (2013) and 
the references therein for their proofs. 

Proposition 2.1 (Prediction). Assume that at a given time t, the pos¬ 
terior cardinality distribution, Pt\t{n), is given and that the posterior PHD 

is a Gaussian mixture of the form Dm{x) = where 

Jt is the number of Gaussian components at t. Then the posterior predicted 
PHD, is also a Gaussian mixture, 

(2-12) Df_|_i|j(x) = bt{x) -\- D, 

where ht{x) is given in (2.11) and Dg t+i\t{x) = PS w^t^N{x-, T%+i\v 

^sl+i\t^ is the PHD which arises from the “survived” organelles. The cor¬ 
responding mean and covariance equal fis,t+i\t = Ftht o-nd Ps,t-i-i\t = Qt + 
FtPtFt", respectively. The posterior predictive cardinality distribution is 

n oo 

(2.13) Pt+i|t(n|Zi,i) = ^pB(n-j)^ 

j=0 l=j 

where pb{-) is the cardinality distribution of the RFS responsible for the 
organelles’ appearance and ps is the survival probability of an organelle. 

We denote the permutations P(f = with the convention that Pff = 

0, if n < m, and we dehne qn = 1 — Pd the probability of not detecting an 
intracellular movement. Furthermore, assume that at time f-|- 1, a new mea¬ 
surement random set, Zt+i, is received with cardinality \Zt+i \ =m. Then 
the predicted PHD (2.12) and cardinality distribution (2.13) will be updated 
according to Proposition 2.2. 

Proposition 2.2 (Update). Suppose that the predicted PHD, Dij_i\f, 
and the cardinality distribution, Pt+iitin^lZi-.t), satisfy Proposition 2.1. Then, 
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the posterior PHD, at a given time t + 1 is a Gaussian mixture, 

and the eorresponding CPHD update equations are listed below: 


- ^Tl+l)- 



Dt+i\tix) 



The mean and the covariance matrix are — 

= V - respectively, where = 

Pt+i\tP^T+i{Pt+i P . Furthermore, the posteriorcardinal- 

ity distribution is propagated via the following equation: 



Remark 2.4. If there were only one intracellular movement during the 
tracking time and neither a birth nor a death of an organelle were allowed, 
then Propositions 2.1 and 2.2 would yield the special case of monitoring a 
random singleton, that is, one organelle in our experiments. Furthermore, if 
the probability of detection pn = I (thus qn = 0) and there was one compo¬ 
nent in the Gaussian mixture, then equation (2.14) would yield the typical 
Kalman filter update equation and in this special case the matrix K would 
play the role of the Kalman gain matrix. 

3. Results. Having established the theoretical framework, we present our 
biological data analysis and tracking in this section. We start with a sum¬ 
mary of our algorithm. 

Step 0: Initialization. The initial intensity, T*o| 0 ) is considered as a Gaus¬ 
sian mixture with Jq components. Furthermore, the initial cardinality dis¬ 
tribution, Po\Q{n), is considered a priori to a single object. 

Step 1: Prediction. At time t the predicted intensity is a Gaussian 

mixture whose components’ weights, means and covariance matrices are de¬ 
rived in equation (2.12). Equation (2.13) yields the corresponding posterior 
predictive cardinality distribution, Pt+i|t(n). 
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Step 2: Update. At time t + 1, the predictions generated in Step 1 are 
updated based on new measurements. More precisely, the posterior PHD, 
is a Gaussian mixture whose weight, mean and covariance ma¬ 
trix is derived by equation (2.14). The posterior cardinality distribution, 
Pt+i\t+i{n), is estimated according to equation (2.15). 

Step 3: Merging and pruning. The number of Gaussian components in¬ 
creases as time progresses. In fact, at a given time, t, the Gaussian mixture 
will require 0{Jt-i\Zt\) components, where Jt-i is the number of compo¬ 
nents of the posterior intensity at time t — 1. Since components with 

low weight do not provide any significant contribution to the approximation 
of the posterior multi-target density, we eliminate the components whose 
weights are negligible and below some preset threshold, T (e.g., T = 10“®). 
The remaining components of the mixture are renormalized such that their 
sum equals 1. 

Furthermore, there are components which are close to each other and prac¬ 
tically could be approximated by a single Gaussian distribution. Indeed, if 
two components of the mixture with weight, state and covariance, {wi, Xi, Pi) 
and {wj,Xj,Pj), respectively, have distance dij = (xj — Xj)P~^{xi — XjY less 
than some threshold, U, then these mixing components are merged into 
one [Glark, Panta and Vo (2006)]. The threshold U should be chosen much 
smaller (e.g., U = 0.004) than the standard deviation of the observations’ 
noise so that the filtering algorithm does not consider two different objects 
as one when they are close together, such as when their paths are crossing 
each other. 

Step 4: Multi-object state extraction. To extract the organelles’ states, 
we focus on only the modes of the corresponding Gaussian mixture. The 
number of organelles is estimated from the cardinality distribution using a 
maximum a posteriori (MAP) estimator h = argsup„p(n|Zi:i). 

Schematically, the algorithm works in the following way, for all t = 0,1,...: 


. ^ ^ Proposition 

where the PHD, D. 
covariance. 


2.1 ,. Proposition 2.2 , x 

|., is estimated via the triplet of weights, mean and 


3.1. Synthetic data. This section illustrates a simulated scenario with 
respect to organelle movements. Consider a set Xt = {x^, ,..., x”‘} whose 

members are 4-dimensional state vectors of the nt organelles at time t. Pre¬ 
cisely, an organelle’s state vector is x\ = \px,t)Vx^t,Py,t,Vy^tY' for any i = 
l,...,nt, where {px,t,Py,t) denote the spatial coordinates of the organelle 
on the xy-plane and the corresponding velocities are denoted as {vx,t-,Vy,t)- 
The movements in a cell may be considered to take place in a force field 
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which is on average inactive. However, when a molecular motor exerts a 
pushing force on an organelle, then there is a positive deviation from the 
mean zero. By the same token, when friction and/or other large enough 
backward-acting forces occur, then the organelles will slow down and even¬ 
tually stop, and thus a symmetric negative deviation from the mean-zero 
force field is caused. Therefore, one may consider that the force field is 
normally distributed with mean zero and pertinent covariance. This con¬ 
sideration is actually validated in Section 3.2 where experimental data are 
analyzed. Given that the mass is conservative over time frames considered 
here (a valid assumption), Newton’s second law yields that the accelera¬ 
tion, a, follows a normal distribution. The velocity changes in turn are also 
normally distributed where the covariance depends on the size of the time 
intervals. Given the fact that Px,t = Vx,t and py^t = Vy,t, we can formally state 
the following linear stochastic differential equations system: 


/dpx,t\ 


/ Vx,t ^ 


/o 



dvx^t 


0 

dt “h 


0 

( dux^t \ 

dPy,t 



0 

0 

V dUy^t ) 

\ dVy^t / 


\ 0 ) 


V 0 

ay) 



where Ux,t,Uy,t are independent Brownian motions and the driving noises 
(JxUx,t and (TyUy^t are Gaussian noises with covariances and ay5{t), 

respectively, where 5{t) is the delta function. Discretizing and approximating 
the system (3.1), we have a two-dimensional model given below: 


(3.2) 


(Px,t ^ 


/I 

A 

0 



iPx,t—l \ 

Vx,t 


0 

1 

0 

0 


Vx,t—l 

Py,t 


0 

0 

1 

A 


Py,t-1 

V Vy^t / 


Vo 

0 

0 

1/ 


V Vy^t-1 / 


/ A2 

T 

A 
0 
V 0 


0 ^ 
0 

A2 


6-1: 


2 

A 


where the model noise, 6-1: is a collection of independent Gaussian random 
variables with covariance matrix S = diag{(T^,Uy}. The sampling time is 
considered A = 1 s since data from organelles’ movements are collected every 
one second. The velocity changes are normally distributed with mean zero, 
and thus 99.7% of the data are within three standard deviations from zero. 
Taking into consideration the biological finding that organelles may move 
up to 7 pm/s (in both directions) [Tominaga et al. (2003)], the standard 
deviation coefficients are chosen ax = (7y = 2.33 pm/s^. If one decreased or 
increased drastically the variance, then the estimates would not be accurate. 
Small noise dynamics (e.g., ax = ay = 0.1 pm/s'^) yield predictions based on 
almost perfect linear dynamics which could lead to erroneous estimation 
in case organelles exhibit a slightly curvy behavior. By the same token, a 
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large standard deviation (e.g., da; = = 5 ^m/s^) produces a wide range of 

samples which lead to inaccurate estimates. 

Furthermore, each object is considered with survival probability, ps,t = 
0.99, such that any organelle within the tracking domain is under monitoring 
unless its signal disappears. The maximum number of involved Gaussian 
components is considered to be fairly large, A^max = 200. The object-birth 
process is a Poisson RFS with intensity defined as in (2.11), where wi, = 0.25, 
= [3 0 5 0]^, /if = [40 -6 0]^,/if) = [-3 0 -2 0]^, /if = [-40 8 0]'^, 

and Pfe = IOI 4 . The four different means, /i^ , i = 1,..., 4 are selected to 
ensure that births on all four quadrants are considered with equal probability 
Wf, = 0.25. The covariance of the birth intensity is also large such that a vast 
candidate area of newborn organelles is covered. Given that our experimental 
environment did not suffer from low signal-to-noise ratio and no false alarm 
occurred, the probability of detecting an organelle is state independent and 
equals pD,t = 0.98. 

We first focus on the synthetic data which consist of the spatial coordi¬ 
nates. Gonsider at given time t -|- 1 the random set, Zt+i = ..., 

where for each i the data = (px,t+i,Py,t+i),i = is a 

two-dimensional vector whose likelihood is defined in ( 2 . 10 ), with 


(3.3) 


Hf = 


1 0 
0 0 


0 0 \ 

1 of 


— 0'1l2 ) 


and (To = 0.2 /rm is the standard deviation of the measurement noise due 
to optical limitations and experimental noise. For example, there is a fun¬ 
damental maximum to the resolution of any optical system due to diffrac¬ 
tion. The diffraction defines the microscope’s point-spread function which 
describes the response of an imaging system to a point light source. Fur¬ 
thermore, our procedure uses a weight threshold T = 10“^ for the pruning 
procedure and a threshold U = 0.004 for the merging part of the algorithm 
(step 3 in the algorithm). 

The synthesized organelles’ trajectories, which play the role of the true 
trajectories, are created by evolving a number of organelles according to 
dynamics (3.2), and the corresponding observations were created after per¬ 
turbing the true trajectories by a normally distributed noise with covariance 
Rt as in (3.3). 

Figure 1 shows that there are twelve organelles (in total) which are mon¬ 
itored for 100 time steps. At any given time t, the number of organelles 
is unknown a priori and is not fixed, that is, random birth and death of 
organelles are allowed with pertinent dynamics based on Assumption 2.3. 
In fact, the organelles’ number increases and decreases drastically during 
the first thirty steps and the last twenty ones as well. This makes the prob¬ 
lem a rather formidable one by keeping in mind that previous studies have 
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Fig. 1. Number of organelles per time step. 


monitored simultaneously a fixed and a priori known number of intracellular 
movements with overall known dynamics, for example, Smal, Niessen and 
Meijering (2006). In contrast, our algorithm assumes an initial cardinality of 
1 (see step 1 of the algorithmic description) and updates its estimate based 
on available data. Thus, our algorithm captures accurately all modifications 
in the number of organelles and it gives an accurate estimate. 

Figure 2 shows a three-dimensional graph of the trajectories’ estimates 
of the organelles across time. As we can see, there are several crossings, 
often in the y-direction. Tracking methods for intracellular movements that 
assume one-to-one correspondence between a measurement and an object 
fail to resolve the most ambiguous track interaction scenarios, for example, 
when objects are in close proximity. However, in our case, we do not assume 



Fig. 2. 


Linear trajectories of organelles in the xy-plane over time. 
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any sort of prior one-to-one correspondence, instead we employ a multi¬ 
object statistical framework by considering a single set of objects, thereby 
producing accurate estimates even during the difficult occasions such as 
crossings. 

Indeed, the estimates are very close to the true trajectories, but to quan¬ 
tify any sort of error a multi-object error distance is considered. The char¬ 
acteristics of a multi-object distance should ( 1 ) be a metric on the space of 
finite sets, (2) capture cardinality and state errors and (3) have a physical 
interpretation. Toward this end, we employ a metric from point processes 
theory in order to measure the discrepancy between the estimates and the 
true values [Bremaud (1981), Mpller and Waagepetersen (2004)]. A formal 
definition of this metric according to Schuhmacher, Vo and Vo (2008) is 
given below. 

Definition 3.1. Let W C be a closed and bounded observation 
window and d denote the Euclidean metric. For c > 0, let d^‘^\x,y) = min(c, 
d{x,y)) denote the distance between x,y and Pn denote the set of 
permutations on {1, 2,..., re} for any re G N. For 1 < i < oo, c > 0 and ar¬ 
bitrary finite subsets X = {xi,... ,Xm} and Y = {yi,... ,y„} of W, where 
m, re = 0,1,2 ,..., define for m<n, 

m \ \ 

and (V, Y) = d!f^ {Y,X) \i m> n. Moreover, if £ = oo, then 
d'^(X,Y) = min max d^^\xi,yT^(i)) if rre = re 

(3.5) 

= c if ire 7 ^ re. 

For any ^ G [1, oo] the distance is equal to zero if m = n = 0. The function 
d!f\x,Y) is called the Optimal SubPattern Assignment (OSPA) metric of 
order i with cutoff parameter c. 

Remark 3.1. Schuhmacher and Xia (2008) examined the special case 
for i = c = 1 and Schuhmacher, Vo and Vo (2008) generalized it for any 

“(c) 

I, c. The metric d} ^ is based on a Wasserstein construction. The advan¬ 
tage of this metric is that equation (3.4) takes into consideration the error 
due to localization and cardinality at the same time. An alternative mea¬ 
sure of discrepancy is the Haussdorff distance [Mpller and Waagepetersen 
(2004)], however, it is relatively insensitive to difference in cardinality as 
was noted in Hoffman and Mahler (2002). The order parameter i is similar 
to the parameter of the Hh order Wasserstein metric between the empirical 


(3,4) rfrD.r)= - mm 
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distributions of the point patterns X and Y. Furthermore, given that c is 
fixed, the parameter i in (3.4) assigns more weight to outliers. The metric 
G [0,c] for any c > 0 in turn gives us a measure of performance 
with respect to the worst possible distance i. Also, if 0 < ci < C 2 < oo, then 
Moreover, the cutoff parameter c determines the weighting of 
how the metric penalizes cardinality errors as opposed to localization errors. 
Smaller values of c tend to put emphasis on the localization error and make 
the metric unchanged by cardinality errors. Thus, the designer can deter¬ 
mine how strongly a false or missing estimate will be penalized by modifying 
the value of c. Here, we have chosen i = \ and c = 30 such that the OSPA 
is sensitive enough in both localization and cardinality errors. The choice of 
the value ^ = 1 has the benefit that the OSPA-metric measures a first order 
per-object error and that the sum of localization and cardinality compo¬ 
nents equals the total metric. The reader should refer to Schuhmacher, Vo 
and Vo (2008) and the references therein for further details on the OSPA 
metric. 


The top picture in Figure 3 depicts the error using the OSPA metric 
given in equation (3.4). We observe that large errors (peaks in the hgure) 
occur when the organelles are crossing and when there is a change of the 
number of organelles (e.g., at t = 20). This is expected since these are the 
most difficult situations. The OSPA error cannot exceed the value 30 since 
the cutoff parameter is set at c = 30, however, even in the most difficult 
cases, the error remains well below 10. The two subsequent pictures are 
showing localization and cardinality error. The localization errors for two 
patterns X = (xi,..., Xm) and Y = (yi,..., y„) with m <n and £ < oo are 
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Fig. 3. Error measured via the OSPA metric. The error can be as large as the cutoff 
parameter c = 30. 
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given by 


4jocD.D= j , 


4%Ax,r) = 


c^{n — m) 


n 


III 


-(c) -(c) 

Strictly speaking, the two errors, and e^^ard’ metrics on the 

space of finite subsets, but one may still gain some insight about the perfor¬ 
mance of the filter [Schuhmacher, Vo and Vo (2008)]. 


3.2. Experimental data. Before outlining our results, we will briefly de¬ 
scribe the conditions under which the movement data were retrieved. Or¬ 
ganelles were labeled with fluorescent protein fusions in root cells of the 
model plant Arabidopsis thaliana and cells on the surface of roots were ob¬ 
served on a fluorescent microscope as described in Nelson, Cai and Nebenfiihr 
(2007). Images were taken with a digital camera at regular intervals (1 s) 
to generate time-lapse sequences of 1 to 2 minute duration (i.e., 60 to 120 
images). These image sequences (e.g.. Figure 4) displayed bright spots of 
different sizes and intensities depending on the size and position of the or¬ 
ganelle relative to the focal plane. Movements of individual organelles were 
readily apparent by comparing the changes in position of spots between im¬ 
age frames (arrow in Figure 4). Specifically, Figure 4 shows the movement of 
peroxisomes, small spherical organelles involved in detoxification of reactive 
oxygen species which have recently emerged as important regulators of plant 
growth and stress responses [Klaus and Heribert (2004)]. Similar movements 
can also be observed for other organelles, such as Golgi stacks [Nebenfiihr 
et al. (1999)] and mitochondria [Van Gestel, Kohler and Verbelen (2002)]. 

Images were analyzed quantitatively by manually marking the center of 
each spot in every frame of the time-lapse sequence which was then recorded 



Fig. 4. Peroxisomes movements. 
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by the Manual Tracking plugin in ImageJ [Schneider, Rasband and Eliceiri 
(2012)]. This procedure produced a series of (x,y) coordinates per image 
frame that were manually linked to specific {x,y) coordinates in subsequent 
frames. The procedure of manually linking is typically slow (about 1 hour 
for the data set that is analyzed herein) and bias due to human decision 
in linking can be a frequent disadvantage. In our case, the resulting two- 
dimensional vectors declaring the position of organelles on the xy-plane at 
every time point were used (1) to calculate the instantaneous velocities of the 
organelles over time; (2) to provide experimental values for the accelerations’ 
distributions; and (3) to provide the raw data to the statistical tracking 
algorithm without knowing a priori which data (coordinates) correspond to 
which organelle. 

In the following, we focus on the motions of eight peroxisomes retrieved 
in experiments in the second author’s lab. First, we decompose the ac¬ 
celeration, and we investigate the distributional behavior of the accelera¬ 
tions per coordinate separately based on the experimental data. There are 
m = 284 acceleration data points from the eight peroxisomes with mean 
and standard deviation on the x-axis, = —0.0326, cj“ = 0.9998, respec¬ 
tively. The corresponding mean and standard deviation on the y-axis are 
y,y = 0.0429, it“ = 0.6922. Next, we test if the accelerations follow a normal 
distribution using a Kolmogorov-Smirnov test and visually by plotting two 
normality plots, one per coordinate. As we can see from the results of the 
Kolmogorov-Smirnov tests presented in Table 1, and the normal probabil¬ 
ity plots in Figure 5, the two accelerations of the eight peroxisomes follow 
a Gaussian distribution. Thus, the arguments of Section 3.1 imply that the 
dynamics of the eight peroxisomes can be described by the discrete system 
in (3.2). 

Therefore, employing the dynamics (3.2) accompanied by the several hy¬ 
perparameters discussed in Section 3.1, we describe our findings for the 
motions of the peroxisomes. Figures 6 and 7 show the trajectories based on 
measurements (line) and the corresponding estimates represented as dots 
in the figures. At the initial time step. Figure 6 shows a greater mismatch 
between the estimates and the data than in the next sampling periods. 

Table 1 

p-values of two Kolmogorov-Smirnov 
tests for the acceleration data points 
of peroxisomes 


Acceleration 

p-value 

Ho 

Oix 

0.31 

Accept 

Oy 

0.3265 

Accept 
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Normal Probability Plot 


Normal Probability Plot 
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Fig. 5. Testing normality of the organelles’ acceleration. Left panel; Acceleration on the 
x-axis. Right panel; Acceleration on the y-axis. 


This is expected since the algorithm attempts to “learn” the pattern of the 
organelles’ motion. Although the peroxisomes’ overall trajectories are not 
linear, they are piecewise linear per time step (1 s), and thus the dynamics 
of Section 3.1 perform satisfactorily since sampling occurs every A = 1 s. 
If the piecewise linearity was violated and/or the acceleration distribution 
was heavy tailed, then the dynamics in equation (3.2) would produce er¬ 
rors which would depend on the curvature of the true trajectories and/or 
the non-Gaussian noise. Figure 8 depicts the cardinality (number of perox¬ 
isomes) per time step. As we observe, the CPHD filter accurately captures 
the target number when their number does not vary, and it takes 1 to 2 
sampling time steps to realize the change in the organelle number. Also, the 
algorithm correctly estimates that there were not any organelles to monitor 
during the time interval [26,29]. The duration of the automated tracking 



Fig. 6. Trajectories of organelles. Left panel; Trajectories in the x-direction. Right panel; 
Trajectories in the y-direction. 




















TRIM: A BAYESIAN APPROACH 


21 





Fig. 7. Trajectories of organelles in the xy-plane over time. 


process based on our algorithm is about 10 s versus roughly 1 hr for the 
manual tracking of the same eight peroxisomes. 

Due to lacking the true trajectories of the organelles (in fact, it is impossi¬ 
ble to know them with the current technology) [Smal, Niessen and Meijering 
(2006)], the OSPA measurement of error (and any other metric of this type) 
cannot be used since it measures the discrepancy between the algorithmic 
estimates and the true trajectories (not the observed measurements). How¬ 
ever, according to our simulation results exposed in Section 3.1, we believe 
that our estimates are very close to the true trajectories of the eight perox¬ 
isomes. 
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Fig. 8. Number of organelles per time step. 
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4. Summary and discussion. In this paper we have considered the motion 
of organelles as evolving sets. This succeeded by incorporating random sets 
techniques for multi-object tracking and using the cardinalized probability 
hypothesis density filter. Employing a novel Gaussian mixture implemen¬ 
tation of the CPHD hlter, we were able to successfully generate an auto¬ 
mated method for a quantitative analysis of intracellular movements, which 
took about 10 seconds versus about 1 hour for manually linking the same 
data. The new approach’s computational cost was linearly dependent on the 
number of objects multiplied by the number of data points. Our model was 
capable of simultaneously monitoring a large number of organelles, specif¬ 
ically peroxisomes, and distinguishing them even when they were in close 
proximity. Consequently, not only did our algorithm monitor the organelles 
but it also gave an accurate estimate on the number of organelles without 
assuming a fixed and known number of them. Furthermore, our data analy¬ 
sis revealed that the acceleration of the peroxisomes are mean-zero normally 
distributed, which according to Newton’s second law supports an on aver¬ 
age “inactive” force field within a cell where positive (pushing force by the 
myosin motors) or backward-acting forces (e.g., friction) are developed in a 
symmetric fashion given that mass is conservative. Consequently, the two 
parameters, myosin power and local friction, were fairly constant on average 
over time and space, respectively. On the other hand, large changes in veloc¬ 
ity (if any) presumably would result from a static organelle engaging with a 
cytoskeletal track, or from a moving organelle dropping from a cytoskeletal 
track. We expect these changes to occur nearly instantaneously, however, 
technical limitations prevented us from detecting these very rapid changes 
if they indeed existed. In particular, we had to employ exposure times up to 
100 ms to obtain sufficient signal for organelle detection. In addition, images 
were taken in 1 s intervals and had a nominal resolution of 200 nm per pixel. 
Given that myosin motors take 35 nm steps and can move up to 7 pm/s, that 
is, one step every 5 ms, as noted in Tominaga et al. (2003), it is apparent 
that these imaging parameters do not allow us to capture the anticipated 
very fast acceleration and deceleration events directly. Instead we can only 
compute the integrated behavior of organelles over many individual myosin 
steps. Therefore, this scientihc conjecture regarding changes in organelle ve¬ 
locities should be further examined on large experimental data sets which 
could yield a more detailed distribution of accelerations, dynamics and thus 
potentially the mechanics within a cell overall. 

Focusing on the algorithm itself, although it captures the organelles’ be¬ 
havior accurately, it did not take other scenarios into consideration which 
would increase the already severe complexity of the problem. For example, 
there might be cases where organelles may move in a more erratic fashion. 
In this scenario, the acceleration distribution might not be normally dis¬ 
tributed and thus nonlinear and/or nonCaussian dynamics could be fruitful 
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for such data. A possible future research avenue is to use high noise with 
suitably controlled drift dynamics or a more complex autoregressive model. 
Another way is to approximate the overall nonlinearities and/or add more 
experimental features, for example, include information about the shape and 
signal intensity of organelles in the linking step [Sbalzarini and Koumout- 
sakos (2005), Smal et al. (2008), Smal, Niessen and Meijering (2006), Smal 
(2009)]. Moreover, the organelles’ survival and detection probabilities were 
presumed state independent and time invariant. On the other hand, these 
probabilities clearly depend on the position of organelles in a cell. For in¬ 
stance, organelles in close proximity to each other may not be detected or, 
given the curvature of cells, the survival probability of an organelle will de¬ 
crease as it approaches an out-of-focus region of the cell. In our experimental 
data, crossings occurred only a few times and organelles were always in-focus 
and “disappeared” when they exited the focal domain. Attempting to by¬ 
pass Assumption 2.2, techniques developed in Hughes, Fricks and Hancock 
(2010), Hughes and Fricks (2011) may be fruitful for these difficult scenarios. 

In conclusion, this manuscript offers the establishment of a systematic 
way of creating an automated algorithm for monitoring motility within a 
cell by considering a unifying statistical framework for multiple objects. 
In turn, such an automated tracking algorithm will greatly strengthen the 
study of motion patterns in cells. Consequently, understanding the typical 
behavior of healthy molecular processes will have a great impact in quickly 
recognizing abnormalities associated with disorders. 
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